##########################################
########################################## 
## STUDY 8
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

#load data for Study 
df <- haven::read_dta(here("Data", "study8.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df_new <- df %>% 
  rename(soc_ladder = q1) %>% 
  rename_at(vars(starts_with("q2_")), ~str_replace(., "q2", "sdrt")) %>%  #status-driven risk-taking
  rename_at(vars(starts_with("q3_")), ~str_replace(., "q3", "dominance")) %>%   #dominance
  rename_at(vars(starts_with("q4_")), ~str_replace(., "q4", "prestige")) %>%   #prestige
  rename_at(vars(starts_with("q5_")), ~str_replace(., "q5", "sdo")) %>%   #social dominance orientation
  rename_at(vars(starts_with("q6_")), ~str_replace(., "q6", "needchaos")) %>%   #need for chaos
  rename_at(vars(starts_with("q7_")), ~str_replace(., "q7", "status")) %>%   #status measures
  rename_at(vars(starts_with("q8_")), ~str_replace(., "q8", "fear_system")) %>%   #fear of opposing system
  rename_at(vars(starts_with("q9_")), ~str_replace(., "q9", "discrmination")) %>%   #discrimination of groups
  rename_at(vars(starts_with("q10")), ~str_replace(., "q10", "sex_orientation")) %>%   #sexual orientation I
  rename_at(vars(starts_with("q11")), ~str_replace(., "q11", "transgender")) %>%   #sexual orientation II
  rename_at(vars(starts_with("q12")), ~str_replace(., "q12", "religion")) %>%   #sexual orientation II
  rename_at(vars(starts_with("q13_")), ~str_replace(., "q13", "att1")) %>%   #attention check I
  rename_at(vars(starts_with("q14_")), ~str_replace(., "q14", "att2")) %>%   #attention check II
  rename_at(vars(starts_with("q15")), ~str_replace(., "q15", "att3")) %>%   #attention check III
  
  #Recode "don't know" as missing
  mutate_at(vars(starts_with("dominance")), ~na_if(., 8)) %>% 
  mutate_at(vars(starts_with("prestige")), ~na_if(., 8)) %>%
  mutate_at(vars(starts_with("needchaos")), ~na_if(., 8)) %>% 
  
  #Construct attention check measure
  mutate(color = tolower(att3),
         color = gsub(" ", "", color),
         color = case_when(
           color == "puce" ~ 1,
           TRUE ~ 0),
         
         black_white = case_when(
           race == 1 ~ "White",
           race == 2 ~ "Black",
           TRUE ~ NA_character_
         ),
         
         male_female = case_when(
           gender == 1 ~ "Women",
           gender == 2 ~ "Men",
           TRUE ~ NA_character_         
           
         ),
         
         #Categorical variable combining race and sex
         race_sex = case_when(
           black_white == "White" & male_female == "Women" ~ "White Women",
           black_white == "White" & male_female == "Men" ~ "White Men",
           black_white == "Black" & male_female == "Men" ~ "Black Men",
           black_white == "Black" & male_female == "Women" ~ "Black Women",
           TRUE ~ NA_character_  
         )) %>% 
  
  #Need for Chaos
  mutate(needchaos_8 = flip(needchaos_8),
         needchaos_9 = flip(needchaos_9)) %>% 
  mutate(nfc = rowMeans(select(., starts_with("needchaos")), na.rm = TRUE)) %>% 
  
  #Status-Driven Risk-Taking
  mutate(sdrt_1 = flip(sdrt_1),
         sdrt_7 = flip(sdrt_7),
         sdrt_8 = flip(sdrt_8),
         sdrt_12 = flip(sdrt_12),
         sdrt_14 = flip(sdrt_14)) %>% 
  mutate(sdrt_scale = rowMeans(select(., starts_with("sdrt_")), na.rm = TRUE)) %>%
  
  #Social Dominance Orientation
  mutate(sdo_3 = flip(sdo_3),
         sdo_4 = flip(sdo_4),
         sdo_7 = flip(sdo_7),
         sdo_8 = flip(sdo_8)) %>% 
  mutate(sdo_scale = rowMeans(select(., starts_with("sdo_")), na.rm = TRUE)) %>% 
  
  #Status measures
  mutate(status_scale = rowMeans(select(., starts_with("status_")), na.rm = TRUE)) %>% 
  
  #Fear of opposing the system
  select(-fear_system_3) %>% 
  mutate(fear_scale = rowMeans(select(., starts_with("fear_system_")), na.rm = TRUE)) %>% 
  
  #Dominance
  mutate(dominance_7 = flip(dominance_7)) %>% 
  mutate(dominance_scale = rowMeans(select(., starts_with("dominance_")), na.rm = TRUE)) %>% 
  
  #Prestige
  mutate(prestige_2 = flip(prestige_2),
         prestige_8 = flip(prestige_8),
         prestige_9 = flip(prestige_9)) %>% 
  mutate(prestige_scale = rowMeans(select(., starts_with("prestige_")), na.rm = TRUE)) %>% 
  
  #Standardize variables
  mutate(across(.cols = c(nfc, sdrt_scale, sdo_scale, status_scale, status_1, status_2, status_3, status_4, fear_scale, dominance_scale, prestige_scale),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_zscore")) %>% 
  
  #Scale variables to range from 0 to 1
  mutate(across(.cols = c(nfc, sdrt_scale, sdo_scale, status_scale, status_1, status_2, status_3, status_4, fear_scale, dominance_scale, prestige_scale),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}_01"))

##########################################
########################################## 
## Include only participants passing color check 
##########################################
########################################## 
df_final <- df_new %>% 
  filter(color == 1)

##########################################
########################################## 
## Sample Descriptives (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
##Demographics among participants identifying as "whites"
df_whites <- df_final %>% 
  filter(black_white == "White")

#gender
print(paste(round(mean(df_whites$gender) - 1, 2)*100 , "percent of participants are Women"))
#age
table(df_whites$profile_age1)
median(df_whites$profile_age1)
#education
print(paste("Less than high school:", round(prop.table(table(df_whites$educ)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df_whites$educ)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df_whites$educ)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df_whites$educ)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df_whites$educ)), 2)[5]*100, "percent"))
print(paste("Graduate or professional degree:", round(prop.table(table(df_whites$educ)), 2)[6]*100, "percent"))

# Demographics among participants identifying as "blacks"
df_blacks<- df_final %>% 
  filter(black_white == "Black")

#gender
print(paste(round(mean(df_blacks$gender) - 1, 2)*100 , "percent of participants are males"))
#age
table(df_blacks$profile_age1)
median(df_blacks$profile_age1)
#education
print(paste("Less than high school:", round(prop.table(table(df_blacks$educ)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df_blacks$educ)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df_blacks$educ)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df_blacks$educ)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df_blacks$educ)), 2)[5]*100, "percent"))
print(paste("Graduate or professional degree:", round(prop.table(table(df_blacks$educ)), 2)[6]*100, "percent"))


##########################################
########################################## 
## Factor Analysis, etc.
##########################################
## Need for Chaos
# Need for Chaos among "Blacks": Cronbach's alpha
nfc_black <- df_final %>% 
  filter(black_white == "Black") %>% 
  select(needchaos_1:needchaos_9) %>% 
  na.omit() %>% 
  psych::alpha(. , check.keys = TRUE) 
nfc_black

print(paste("Cronbach's alpha, Need for Chaos:" , round(nfc_black$total$raw_alpha,2)))  

# Need for Chaos among "Whites": Cronbach's alpha
nfc_white <- df_final %>% 
  filter(black_white == "White") %>% 
  select(needchaos_1:needchaos_9) %>% 
  na.omit() %>% 
  psych::alpha(. , check.keys = TRUE) 
nfc_white
print(paste("Cronbach's alpha, Need for Chaos:" , round(nfc_white$total$raw_alpha,2)))  

#Confirmatory factor analysis, Need for Chaos: All participants
cfa <- df_final %>% 
  select(needchaos_1:needchaos_9) %>% 
  na.omit()

#One-factor model (Table S7)
m1  <- ' nfc_oneDim  =~ needchaos_1 + needchaos_2 + needchaos_3 + needchaos_4 + needchaos_5 + needchaos_6 + needchaos_7 + needchaos_8 + needchaos_9
        needchaos_4~~needchaos_5
        needchaos_4~~needchaos_6
        needchaos_5~~needchaos_6
        needchaos_8~~needchaos_9'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

## Status-Driven Risk-Taking: Cronbach's alpha
sdrt <- df_final %>% 
  select(sdrt_1:sdrt_14) %>% 
  psych::alpha(. , check.keys = TRUE)
sdrt

## Social Dominance Orientation: Cronbach's alpha
sdo <- df_final %>% 
  select(sdo_1:sdo_8) %>% 
  psych::alpha(. , check.keys = TRUE)
sdo

## Status Motivations: Cronbach's alpha
status <- df_final %>% 
  select(status_1:status_4) %>% 
  psych::alpha(. , check.keys = TRUE)
status

## Fear of Opposing the System: Cronbach's alpha
fear <- df_final %>% 
  select(fear_system_1:fear_system_4) %>% 
  psych::alpha(. , check.keys = TRUE)
fear

## Dominance: Cronbach's alpha
dom <- df_final %>% 
  select(dominance_1:dominance_8) %>% 
  psych::alpha(. , check.keys = TRUE)
dom

## Prestige: Cronbach's alpha
pres <- df_final %>% 
  select(prestige_1:prestige_9) %>% 
  psych::alpha(. , check.keys = TRUE)
pres

##########################################
########################################## 
## Remove 5% with highest Need for Chaos? (SM8)
## When set to TRUE, used for producing Fig S13, Fig S14, Fig S15
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df_final <- df_final %>% 
    filter(nfc_01 < quantile(nfc_01, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 8
##########################################
########################################## 
df_reduced <- df_final %>% 
  select(status_1_01:status_4_01, fear_scale_01,  sdrt_scale_01, nfc_01, dominance_scale_01, prestige_scale_01,  profile_age1,  
         educ, race_sex, faminc_new, black_white, male_female, fear_scale_01, status_scale_01, status_1_zscore:status_4_zscore ) %>% 
  na.omit() %>% data.frame()

#Regression models of Need for Chaos on four status measures in same model and separately
reg1 <- lm(nfc_01 ~ status_1_zscore + status_2_zscore + status_3_zscore + status_4_zscore + factor(educ) + profile_age1 + factor(race_sex) + factor(faminc_new) , 
           data = df_reduced)
reg2 <- lm(nfc_01 ~ status_1_zscore + factor(educ) + profile_age1 + factor(race_sex) + factor(faminc_new) , 
           data = df_reduced)
reg3 <- lm(nfc_01 ~ status_2_zscore + factor(educ) + profile_age1 + factor(race_sex) + factor(faminc_new) , 
           data = df_reduced)
reg4 <- lm(nfc_01 ~ status_3_zscore + factor(educ) + profile_age1 + factor(race_sex) + factor(faminc_new) , 
           data = df_reduced)
reg5 <- lm(nfc_01 ~ status_4_zscore + factor(educ) + profile_age1 + factor(race_sex) + factor(faminc_new) , 
           data = df_reduced)

# Plot estimates from aggregate model
df_plot <- tidy(reg1) %>% 
  filter(term %in% c("status_1_zscore", "status_2_zscore", "status_3_zscore", "status_4_zscore"))

df_plot$term <- factor(df_plot$term, levels = df_plot$term[order(df_plot$estimate)])
df_plot$`Status Type` <- plyr::revalue(df_plot$term , c("status_1_zscore" = "Personal Status: Gain", 
                                                        "status_2_zscore" = "Group Status: Gain", 
                                                        "status_3_zscore" = "Personal Status: Loss",
                                                        "status_4_zscore" = "Group Status: Loss"))

df_plot <- df_plot %>% 
  mutate(`Status Type` = factor(`Status Type`)) %>% 
  mutate(`Status Type`  = fct_relevel(`Status Type`, "Group Status: Loss", "Group Status: Gain", 
                                      "Personal Status: Gain","Personal Status: Loss"))

text <- paste("beta == ", round(df_plot$estimate , digits = 2) )

fig7 <- ggplot(df_plot) +
  geom_hline(yintercept = 0, color = "black") +
  geom_point(aes(x = term , y = estimate, shape = `Status Type`), size = 3) +
  geom_linerange(aes(x = term, y = estimate, ymin = estimate-2*std.error, ymax = estimate+2*std.error), 
                 size = 2, alpha = .2) +
  coord_flip() + 
  theme_bw() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.title.x = element_text(size = 16),
    plot.title = element_text(face="bold")) +
  labs(title = ""  ) +
  ylim(-.1,.2) +
  xlab("") +
  ylab("") +
  theme(legend.justification=c(1,0), legend.position=c(.98,0.05), legend.text=element_text(size=12)) +
  ylab("Estimated Regression Coefficients: Status on Need for Chaos") +
  geom_text(aes(x = term, y = estimate),
            label=text, 
            nudge_y = 0.0, nudge_x = .4, 
            check_overlap = T , parse = TRUE, size = 6
  ) +
  guides(fill = guide_legend(reverse=TRUE), shape = guide_legend(reverse = TRUE)) +
  scale_y_continuous(breaks = c(-.1,-.05,0,.05,.1,.15,.2,.25,.3), 
                     labels = c("-.1","","0","",".1","",".2","",".3"),
                     limits = c(-.1,.2))

ggsave(here::here("figures", "fig8.pdf"),
       width = 7, height = 6)

# Regression Table for Supplemental Materials (Table S22)
stargazer::stargazer(reg1,reg2,reg3,reg4,reg5,
                     align=TRUE,
                     no.space=TRUE,
                     covariate.labels = c("Personal Status: Gain (std.)", "Group Status: Gain (std.)" , "Personal Status: Loss (std.)" , 
                                          "Group Status: Loss (std.)", "High school graduate (ref: No HS)" ,"Some College", "2-year",
                                          "4-year", "Post-Grad", "Age",
                                          "Black Women (ref: Black Men)", "White Men (ref: Black Men)",
                                          "White Women Women (ref: Black Men)", 
                                          "HH Income: $10,000 - $19,999 (ref: < $10,000)", "HH Income: $20,000 - $29,999",
                                          "HH Income: $30,000 - $39,999","HH Income: $40,000 - $49,999","HH Income: $50,000 - $59,999",
                                          "HH Income: $60,000 - $69,999","HH Income: $70,000 - $79,999",
                                          "HH Income: $80,000 - $99,999","HH Income: $100,000 - $119,999",
                                          "HH Income: $120,000 - $149,999","HH Income: $150,000 - $199,999",
                                          "HH Income: $200,000 - $249,999","HH Income: $250,000 - $349,999",
                                          "HH Income: $350,000 - $499,999","HH Income: $500,000 or more",
                                          "HH Income: Prefer not to say","Intercept"),
                     keep.stat = c("n"))

##########################################
########################################## 
## Analysis - Fig. 9
##########################################
########################################## 
## Panel A of Fig. 9: Need for Chaos
a <- plot_fun(df = df_final,
              x = race_sex,
              y = nfc_01,
              text = "does race/gender predict need for chaos?")
nfc_fig8 <- rbind(a[[1]])

fig8_a <- ggplot(nfc_fig8) +
  geom_point(aes(x = race_sex, y = mean), size = 3,
             position = position_dodge(.55)) +
  geom_linerange(aes(x = race_sex,y = mean, ymax = high_ci, 
                     ymin = low_ci,),
                 position = position_dodge(.55),
                 size =2.5,
                 alpha = .2) +
  ylab("Need for Chaos") + xlab("") +
  ylim(c(0.15,.45)) +
  coord_flip() + 
  theme_bw() +
  theme(
    legend.position = "none",
    axis.text.x = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.text.y = element_text(size = 13),
    plot.title = element_text(face="bold")) +
  ggtitle("A: Need for Chaos")

## Panel B of Fig. 9: Status Measures
b <- plot_fun(df = df_final,
              x = race_sex,
              y = status_1_01,
              text = "does race/gender predict personal status gain?")

c <- plot_fun(df = df_final,
              x = race_sex,
              y = status_2_01,
              text = "does race/gender predict group status gain?")

d <- plot_fun(df = df_final,
              x = race_sex,
              y = status_3_01,
              text = "does race/gender predict personal statuss loss?")

e <- plot_fun(df = df_final,
              x = race_sex,
              y = status_4_01,
              text = "does race/gender predict group status loss?")

status_fig8 <- rbind(b[[1]], c[[1]], d[[1]], e[[1]])

status_fig8 <- status_fig8 %>% 
  mutate(`Status Type` = rep(c("Personal Status Gain", "Group Status Gain", "Personal Status Loss", "Group Status Loss"), each = 4),
         race_sex = factor(race_sex)) %>% 
  tibble(.)

#Dataframe for labels
fig_labels <-  status_fig8 %>%
  filter(row_number() %in% c(4,8,12,16))

fig8_b <- ggplot(status_fig8) +
  geom_point(aes(x = race_sex, y = mean, fill = `Status Type`, shape = `Status Type`), size = 3,
             position = position_dodge(.55)) +
  geom_linerange(aes(x = race_sex,y = mean, ymax = high_ci, 
                     ymin = low_ci, fill = `Status Type`),
                 position = position_dodge(.55),
                 size =2.5,
                 alpha = .2) +
  ylab("Status Concerns") + xlab("") +
  ylim(c(0.2,.8)) +
  coord_flip() + 
  theme_bw() +
  theme(
    legend.position=c(.72,.8),
    legend.text=element_text(size=12),
    axis.text.x = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.text.y = element_text(size = 13),
    plot.title = element_text(face="bold")) +
  scale_x_discrete(position = "top") +
  ggtitle("B: Status Concerns") +
  guides(fill = guide_legend(reverse=TRUE), shape = guide_legend(reverse = TRUE))

#Add plots
figure <- ggarrange(fig8_a, fig8_b,
                    ncol = 2, nrow = 1)
#Save plot
ggsave(here::here("Figures", "fig9.pdf"),
       width = 10.5, height =5.5)

## Regression Table for Supplemental Materials (Table S23 & S24)
m1 <- lm(nfc_01 ~ -1 + factor(male_female), data = df_final)
m2 <- lm(nfc_01 ~ -1 + factor(black_white), data = df_final)
m3 <- lm(nfc_01 ~ -1 + factor(race_sex), data = df_final)

stargazer::stargazer(m1,m2,m3,
                     align=TRUE,
                     no.space=TRUE,
                     covariate.labels = c("Women", "Men", "Black", "White", "Black Men", "Black Women", "White Men", "White Women" ,"Intercept"),
                     keep.stat = c("n"), ci = TRUE, ci.level = 0.95, star.cutoffs = NA, omit.table.layout = "n",
                     p = NULL)

m1 <- lm(status_1_01 ~ -1 + factor(race_sex), data = df_final)
m2 <- lm(status_2_01 ~ -1 + factor(race_sex), data = df_final)
m3 <- lm(status_3_01 ~ -1 + factor(race_sex), data = df_final)
m4 <- lm(status_4_01 ~ -1 + factor(race_sex), data = df_final)

stargazer::stargazer(m1,m2,m3,m4,
                     align=TRUE,
                     no.space=TRUE,
                     dep.var.labels = c("Personal Gain", "Group Gain", "Personal Loss", "Group Loss"),
                     covariate.labels = c("Black Men", "Black Women", "White Men", "White Women"),
                     keep.stat = c("n"), ci = TRUE, ci.level = 0.95, star.cutoffs = NA, omit.table.layout = "n",
                     p = NULL)

##########################################
########################################## 
## Analysis - Fig. 10
##########################################
########################################## 
## Status 1
m1 <- lm(nfc_01 ~ status_1_zscore*race_sex  + status_2_zscore + status_3_zscore + status_4_zscore, data = df_final)
m <- margins::margins(m1)
spl <- split(m, m$race_sex)

white_women_df1 <- tidy(spl$`White Women`) %>% 
  filter(term %in% c("status_1_zscore"))
white_men_df1 <- tidy(spl$`White Men`)  %>% 
  filter(term %in% c("status_1_zscore"))
black_women_df1 <- tidy(spl$`Black Women`)  %>% 
  filter(term %in% c("status_1_zscore"))
black_men_df1 <- tidy(spl$`Black Men`)  %>% 
  filter(term %in% c("status_1_zscore"))
## Status 2
m2 <- lm(nfc_01 ~ status_2_zscore*race_sex  + status_3_zscore + status_1_zscore + status_4_zscore, data = df_final)
m <- margins::margins(m2)
spl <- split(m, m$race_sex)

white_women_df2 <- tidy(spl$`White Women`)  %>% 
  filter(term %in% c("status_2_zscore"))
white_men_df2 <- tidy(spl$`White Men`)  %>% 
  filter(term %in% c("status_2_zscore"))
black_women_df2 <- tidy(spl$`Black Women`)  %>% 
  filter(term %in% c("status_2_zscore"))
black_men_df2 <- tidy(spl$`Black Men`)  %>% 
  filter(term %in% c("status_2_zscore"))
## Status 3
m3 <- lm(nfc_01 ~ status_3_zscore*race_sex + status_2_zscore + status_1_zscore + status_4_zscore, data = df_final)
m <- margins::margins(m3)
spl <- split(m, m$race_sex)

white_women_df3 <- tidy(spl$`White Women`)  %>% 
  filter(term %in% c("status_3_zscore"))
white_men_df3 <- tidy(spl$`White Men`)  %>% 
  filter(term %in% c("status_3_zscore"))
black_women_df3 <- tidy(spl$`Black Women`)  %>% 
  filter(term %in% c("status_3_zscore"))
black_men_df3 <- tidy(spl$`Black Men`)  %>% 
  filter(term %in% c("status_3_zscore"))
## Status 4
m4 <- lm(nfc_01 ~ status_4_zscore*race_sex  + status_2_zscore + status_1_zscore + status_3_zscore, data = df_final)
m <- margins::margins(m4)
spl <- split(m, m$race_sex)

white_women_df4 <- tidy(spl$`White Women`)  %>% 
  filter(term %in% c("status_4_zscore"))
white_men_df4 <- tidy(spl$`White Men`)  %>% 
  filter(term %in% c("status_4_zscore"))
black_women_df4 <- tidy(spl$`Black Women`)  %>% 
  filter(term %in% c("status_4_zscore"))
black_men_df4 <- tidy(spl$`Black Men`)  %>% 
  filter(term %in% c("status_4_zscore"))

df_all <- rbind(white_women_df1,white_women_df2,white_women_df3,white_women_df4,
                white_men_df1,white_men_df2,white_men_df3,white_men_df4,
                black_women_df1,black_women_df2,black_women_df3,black_women_df4,
                black_men_df1,black_men_df2,black_men_df3,black_men_df4) %>% 
  mutate(`Demographic Group` = rep(c("White Women", "White Men", "Black Women", "Black Men"), each = 4),
         `Demographic Group` = factor(`Demographic Group`),
         `Demographic Group` = fct_relevel(`Demographic Group`, "Black Women", "Black Men", "White Women","White Men" )) %>% 
  mutate(`Status Type` = rep(c("Personal Status Gain", "Group Status Gain", "Personal Status Loss", "Group Status Loss"), times = 4))

fig9 <- ggplot(df_all) +
  geom_hline(yintercept = 0, color = "black") +
  geom_point(aes(x = `Status Type`, y = estimate, fill = `Demographic Group`, shape = `Demographic Group`), size = 3,
             position = position_dodge(.55)) +
  geom_linerange(aes(x = `Status Type`,y = estimate, ymax = estimate + 2*std.error, 
                     ymin = estimate - 2*std.error, fill = `Demographic Group`),
                 position = position_dodge(.55),
                 size =2.5,
                 alpha = .2) +
  ylab("Association Between Status Concerns and Need for Chaos") + xlab("") +
  ylim(c(-0.3,.5)) +
  coord_flip() + 
  theme_bw() +
  theme(
    legend.position=c(.8,.2),
    axis.text.x = element_text(size = 13),
    legend.text=element_text(size=12),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    plot.title = element_text(face="bold", size = 12)) +
  ggtitle("") +
  guides(fill = guide_legend(reverse=TRUE), shape = guide_legend(reverse = TRUE)) +
  scale_y_continuous(n.breaks = 8, limits = c(-.1,.2))

ggsave(here::here("Figures", "fig10.pdf"),
       width = 8, height =5.5)

## Regression Table for Supplemental Materials (Table S25)
m1 <- lm(nfc_01 ~ status_1_zscore*relevel(as.factor(race_sex), ref = "White Men") + status_2_zscore + status_3_zscore + status_4_zscore, data = df_final)
m2 <- lm(nfc_01 ~ status_2_zscore*relevel(as.factor(race_sex), ref = "White Men") + status_1_zscore + status_3_zscore + status_4_zscore, data = df_final)
m3 <- lm(nfc_01 ~ status_3_zscore*relevel(as.factor(race_sex), ref = "White Men") + status_2_zscore + status_1_zscore + status_4_zscore, data = df_final)
m4 <- lm(nfc_01 ~ status_4_zscore*relevel(as.factor(race_sex), ref = "White Men") + status_2_zscore + status_3_zscore + status_1_zscore, data = df_final)

stargazer(m1,m2,m3, m4,
          align=TRUE,
          no.space=TRUE,
          keep.stat = c("n"),
          covariate.labels = c("Personal Gain (PG)", "GL X Black Men", "GL X Black Women", "GL X White Women", 
                               "Black Men (ref. White Men)", "Black Women (ref. White Men)",
                               "White Women (ref. White Men)", "Group Gain (GG)", "Personal Loss (PL)", "Group Loss (GL)", 
                               "PG X Black Men", "PG X Black Women", "PG X White Women", 
                               "GG X Black Men", "GG X Black Women", "GG X White Women",
                               "PL X Black Men", "PL X Black Women", "PL X White Women",  
                               "Intercept"
          ),
          dep.var.labels.include = FALSE,
          font.size = "small", 
          omit.table.layout = "n")

##########################################
########################################## 
## Analysis - Fig S5
##########################################
########################################## 
s1.binning <- interflex(Y = "nfc_01", X = "status_1_zscore", D="race_sex", Z = c("status_2_zscore", "status_3_zscore", "status_4_zscore"),
                        data = df_reduced,  
                        estimator = "linear")
a <- predict(s1.binning, xlab = "Personal Status: Gain", ylab = "Need for Chaos") 

s2.binning <- interflex(Y = "nfc_01", X = "status_2_zscore", D="race_sex", Z = c("status_1_zscore", "status_3_zscore", "status_4_zscore"),
                        data = df_reduced,  
                        estimator = "linear")
b <- predict(s2.binning, xlab = "Group Status: Gain", ylab = "Need for Chaos") 

s3.binning <- interflex(Y = "nfc_01", X = "status_3_zscore", D="race_sex", Z = c("status_2_zscore", "status_1_zscore", "status_4_zscore"),
                        data = df_reduced,  
                        estimator = "linear")
c <- predict(s3.binning, xlab = "Personal Status: Loss", ylab = "Need for Chaos") 


s4.binning <- interflex(Y = "nfc_01", X = "status_4_zscore", D="race_sex", Z = c("status_2_zscore", "status_3_zscore", "status_1_zscore"),
                        data = df_reduced,  
                        estimator = "linear")
d <- predict(s4.binning, xlab = "Group Status: Loss", ylab = "Need for Chaos") 

ggarrange(a,b,c,d,
          ncol = 1, nrow = 4)

ggsave(here::here("figures", "figS5.pdf"),
       width = 10, height =12)

##########################################
########################################## 
## Analysis - Fig S1
##########################################
########################################## 
df_corr <- df_final %>% 
  select(nfc_01, prestige_scale_01, dominance_scale_01, sdrt_scale_01) %>% 
  rename(`Need for Chaos` = nfc_01,
         `Prestige` = prestige_scale_01,
         `Dominance` = dominance_scale_01,
         `Status-Driven Risk-Taking` = sdrt_scale_01) %>% 
  na.omit()


pdf(file = here::here("figures", "figS1.pdf"),   
    width = 18,
    height = 12) 

chart.Correlation(df_corr, histogram=TRUE, cex = 5)

dev.off()

##########################################
########################################## 
## Analysis - Fig S6
##########################################
########################################## 
a <- plot_fun(df = df_final,
              x = race_sex,
              y = fear_scale_01,
              text = "Sociodemographic groups and fear of opposing system")

a[[2]] + ylab("Fear of Opposing System")

ggsave(here::here("figures", "figS6.pdf"),
       width = 8, height =6)

## Regression Table for Supplemental Materials (Table S26)
m <- lm(nfc_01 ~ fear_scale_01*relevel(factor(race_sex), ref = "Black Men") + fear_scale_01*educ + fear_scale_01*profile_age1 + 
          fear_scale_01*faminc_new,
        data = df_reduced)

stargazer(m,
          align=TRUE,
          no.space=TRUE,
          keep.stat = c("n"),
          covariate.labels = c("Fear of Opposing System (0-1)" , "Black Women (ref.: Black Men)", "White Women (ref.: Black Men)", "White Men (ref.: Black Men)", "Education",
                               "Age", "Income", "Fear of System X Black Women", "Fear of System X White Men", "Fear of System X White Women",
                               "Fear of System X Education", "Fear of System X Age", "Fear of System X Income", "Intercept"),
          dep.var.labels.include = FALSE)







